home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / zlabrd.f < prev    next >
Text File  |  1996-07-19  |  13KB  |  330 lines
  1.       SUBROUTINE ZLABRD( M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y,
  2.      $                   LDY )
  3. *
  4. *  -- LAPACK auxiliary routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       INTEGER            LDA, LDX, LDY, M, N, NB
  11. *     ..
  12. *     .. Array Arguments ..
  13.       DOUBLE PRECISION   D( * ), E( * )
  14.       COMPLEX*16         A( LDA, * ), TAUP( * ), TAUQ( * ), X( LDX, * ),
  15.      $                   Y( LDY, * )
  16. *     ..
  17. *
  18. *  Purpose
  19. *  =======
  20. *
  21. *  ZLABRD reduces the first NB rows and columns of a complex general
  22. *  m by n matrix A to upper or lower real bidiagonal form by a unitary
  23. *  transformation Q' * A * P, and returns the matrices X and Y which
  24. *  are needed to apply the transformation to the unreduced part of A.
  25. *
  26. *  If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower
  27. *  bidiagonal form.
  28. *
  29. *  This is an auxiliary routine called by ZGEBRD
  30. *
  31. *  Arguments
  32. *  =========
  33. *
  34. *  M       (input) INTEGER
  35. *          The number of rows in the matrix A.
  36. *
  37. *  N       (input) INTEGER
  38. *          The number of columns in the matrix A.
  39. *
  40. *  NB      (input) INTEGER
  41. *          The number of leading rows and columns of A to be reduced.
  42. *
  43. *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
  44. *          On entry, the m by n general matrix to be reduced.
  45. *          On exit, the first NB rows and columns of the matrix are
  46. *          overwritten; the rest of the array is unchanged.
  47. *          If m >= n, elements on and below the diagonal in the first NB
  48. *            columns, with the array TAUQ, represent the unitary
  49. *            matrix Q as a product of elementary reflectors; and
  50. *            elements above the diagonal in the first NB rows, with the
  51. *            array TAUP, represent the unitary matrix P as a product
  52. *            of elementary reflectors.
  53. *          If m < n, elements below the diagonal in the first NB
  54. *            columns, with the array TAUQ, represent the unitary
  55. *            matrix Q as a product of elementary reflectors, and
  56. *            elements on and above the diagonal in the first NB rows,
  57. *            with the array TAUP, represent the unitary matrix P as
  58. *            a product of elementary reflectors.
  59. *          See Further Details.
  60. *
  61. *  LDA     (input) INTEGER
  62. *          The leading dimension of the array A.  LDA >= max(1,M).
  63. *
  64. *  D       (output) DOUBLE PRECISION array, dimension (NB)
  65. *          The diagonal elements of the first NB rows and columns of
  66. *          the reduced matrix.  D(i) = A(i,i).
  67. *
  68. *  E       (output) DOUBLE PRECISION array, dimension (NB)
  69. *          The off-diagonal elements of the first NB rows and columns of
  70. *          the reduced matrix.
  71. *
  72. *  TAUQ    (output) COMPLEX*16 array dimension (NB)
  73. *          The scalar factors of the elementary reflectors which
  74. *          represent the unitary matrix Q. See Further Details.
  75. *
  76. *  TAUP    (output) COMPLEX*16 array, dimension (NB)
  77. *          The scalar factors of the elementary reflectors which
  78. *          represent the unitary matrix P. See Further Details.
  79. *
  80. *  X       (output) COMPLEX*16 array, dimension (LDX,NB)
  81. *          The m-by-nb matrix X required to update the unreduced part
  82. *          of A.
  83. *
  84. *  LDX     (input) INTEGER
  85. *          The leading dimension of the array X. LDX >= max(1,M).
  86. *
  87. *  Y       (output) COMPLEX*16 array, dimension (LDY,NB)
  88. *          The n-by-nb matrix Y required to update the unreduced part
  89. *          of A.
  90. *
  91. *  LDY     (output) INTEGER
  92. *          The leading dimension of the array Y. LDY >= max(1,N).
  93. *
  94. *  Further Details
  95. *  ===============
  96. *
  97. *  The matrices Q and P are represented as products of elementary
  98. *  reflectors:
  99. *
  100. *     Q = H(1) H(2) . . . H(nb)  and  P = G(1) G(2) . . . G(nb)
  101. *
  102. *  Each H(i) and G(i) has the form:
  103. *
  104. *     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u'
  105. *
  106. *  where tauq and taup are complex scalars, and v and u are complex
  107. *  vectors.
  108. *
  109. *  If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in
  110. *  A(i:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is stored on exit in
  111. *  A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
  112. *
  113. *  If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in
  114. *  A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i:n) is stored on exit in
  115. *  A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
  116. *
  117. *  The elements of the vectors v and u together form the m-by-nb matrix
  118. *  V and the nb-by-n matrix U' which are needed, with X and Y, to apply
  119. *  the transformation to the unreduced part of the matrix, using a block
  120. *  update of the form:  A := A - V*Y' - X*U'.
  121. *
  122. *  The contents of A on exit are illustrated by the following examples
  123. *  with nb = 2:
  124. *
  125. *  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):
  126. *
  127. *    (  1   1   u1  u1  u1 )           (  1   u1  u1  u1  u1  u1 )
  128. *    (  v1  1   1   u2  u2 )           (  1   1   u2  u2  u2  u2 )
  129. *    (  v1  v2  a   a   a  )           (  v1  1   a   a   a   a  )
  130. *    (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  )
  131. *    (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  )
  132. *    (  v1  v2  a   a   a  )
  133. *
  134. *  where a denotes an element of the original matrix which is unchanged,
  135. *  vi denotes an element of the vector defining H(i), and ui an element
  136. *  of the vector defining G(i).
  137. *
  138. *  =====================================================================
  139. *
  140. *     .. Parameters ..
  141.       COMPLEX*16         ZERO, ONE
  142.       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ),
  143.      $                   ONE = ( 1.0D+0, 0.0D+0 ) )
  144. *     ..
  145. *     .. Local Scalars ..
  146.       INTEGER            I
  147.       COMPLEX*16         ALPHA
  148. *     ..
  149. *     .. External Subroutines ..
  150.       EXTERNAL           ZGEMV, ZLACGV, ZLARFG, ZSCAL
  151. *     ..
  152. *     .. Intrinsic Functions ..
  153.       INTRINSIC          MIN
  154. *     ..
  155. *     .. Executable Statements ..
  156. *
  157. *     Quick return if possible
  158. *
  159.       IF( M.LE.0 .OR. N.LE.0 )
  160.      $   RETURN
  161. *
  162.       IF( M.GE.N ) THEN
  163. *
  164. *        Reduce to upper bidiagonal form
  165. *
  166.          DO 10 I = 1, NB
  167. *
  168. *           Update A(i:m,i)
  169. *
  170.             CALL ZLACGV( I-1, Y( I, 1 ), LDY )
  171.             CALL ZGEMV( 'No transpose', M-I+1, I-1, -ONE, A( I, 1 ),
  172.      $                  LDA, Y( I, 1 ), LDY, ONE, A( I, I ), 1 )
  173.             CALL ZLACGV( I-1, Y( I, 1 ), LDY )
  174.             CALL ZGEMV( 'No transpose', M-I+1, I-1, -ONE, X( I, 1 ),
  175.      $                  LDX, A( 1, I ), 1, ONE, A( I, I ), 1 )
  176. *
  177. *           Generate reflection Q(i) to annihilate A(i+1:m,i)
  178. *
  179.             ALPHA = A( I, I )
  180.             CALL ZLARFG( M-I+1, ALPHA, A( MIN( I+1, M ), I ), 1,
  181.      $                   TAUQ( I ) )
  182.             D( I ) = ALPHA
  183.             IF( I.LT.N ) THEN
  184.                A( I, I ) = ONE
  185. *
  186. *              Compute Y(i+1:n,i)
  187. *
  188.                CALL ZGEMV( 'Conjugate transpose', M-I+1, N-I, ONE,
  189.      $                     A( I, I+1 ), LDA, A( I, I ), 1, ZERO,
  190.      $                     Y( I+1, I ), 1 )
  191.                CALL ZGEMV( 'Conjugate transpose', M-I+1, I-1, ONE,
  192.      $                     A( I, 1 ), LDA, A( I, I ), 1, ZERO,
  193.      $                     Y( 1, I ), 1 )
  194.                CALL ZGEMV( 'No transpose', N-I, I-1, -ONE, Y( I+1, 1 ),
  195.      $                     LDY, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
  196.                CALL ZGEMV( 'Conjugate transpose', M-I+1, I-1, ONE,
  197.      $                     X( I, 1 ), LDX, A( I, I ), 1, ZERO,
  198.      $                     Y( 1, I ), 1 )
  199.                CALL ZGEMV( 'Conjugate transpose', I-1, N-I, -ONE,
  200.      $                     A( 1, I+1 ), LDA, Y( 1, I ), 1, ONE,
  201.      $                     Y( I+1, I ), 1 )
  202.                CALL ZSCAL( N-I, TAUQ( I ), Y( I+1, I ), 1 )
  203. *
  204. *              Update A(i,i+1:n)
  205. *
  206.                CALL ZLACGV( N-I, A( I, I+1 ), LDA )
  207.                CALL ZLACGV( I, A( I, 1 ), LDA )
  208.                CALL ZGEMV( 'No transpose', N-I, I, -ONE, Y( I+1, 1 ),
  209.      $                     LDY, A( I, 1 ), LDA, ONE, A( I, I+1 ), LDA )
  210.                CALL ZLACGV( I, A( I, 1 ), LDA )
  211.                CALL ZLACGV( I-1, X( I, 1 ), LDX )
  212.                CALL ZGEMV( 'Conjugate transpose', I-1, N-I, -ONE,
  213.      $                     A( 1, I+1 ), LDA, X( I, 1 ), LDX, ONE,
  214.      $                     A( I, I+1 ), LDA )
  215.                CALL ZLACGV( I-1, X( I, 1 ), LDX )
  216. *
  217. *              Generate reflection P(i) to annihilate A(i,i+2:n)
  218. *
  219.                ALPHA = A( I, I+1 )
  220.                CALL ZLARFG( N-I, ALPHA, A( I, MIN( I+2, N ) ), LDA,
  221.      $                      TAUP( I ) )
  222.                E( I ) = ALPHA
  223.                A( I, I+1 ) = ONE
  224. *
  225. *              Compute X(i+1:m,i)
  226. *
  227.                CALL ZGEMV( 'No transpose', M-I, N-I, ONE, A( I+1, I+1 ),
  228.      $                     LDA, A( I, I+1 ), LDA, ZERO, X( I+1, I ), 1 )
  229.                CALL ZGEMV( 'Conjugate transpose', N-I, I, ONE,
  230.      $                     Y( I+1, 1 ), LDY, A( I, I+1 ), LDA, ZERO,
  231.      $                     X( 1, I ), 1 )
  232.                CALL ZGEMV( 'No transpose', M-I, I, -ONE, A( I+1, 1 ),
  233.      $                     LDA, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
  234.                CALL ZGEMV( 'No transpose', I-1, N-I, ONE, A( 1, I+1 ),
  235.      $                     LDA, A( I, I+1 ), LDA, ZERO, X( 1, I ), 1 )
  236.                CALL ZGEMV( 'No transpose', M-I, I-1, -ONE, X( I+1, 1 ),
  237.      $                     LDX, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
  238.                CALL ZSCAL( M-I, TAUP( I ), X( I+1, I ), 1 )
  239.                CALL ZLACGV( N-I, A( I, I+1 ), LDA )
  240.             END IF
  241.    10    CONTINUE
  242.       ELSE
  243. *
  244. *        Reduce to lower bidiagonal form
  245. *
  246.          DO 20 I = 1, NB
  247. *
  248. *           Update A(i,i:n)
  249. *
  250.             CALL ZLACGV( N-I+1, A( I, I ), LDA )
  251.             CALL ZLACGV( I-1, A( I, 1 ), LDA )
  252.             CALL ZGEMV( 'No transpose', N-I+1, I-1, -ONE, Y( I, 1 ),
  253.      $                  LDY, A( I, 1 ), LDA, ONE, A( I, I ), LDA )
  254.             CALL ZLACGV( I-1, A( I, 1 ), LDA )
  255.             CALL ZLACGV( I-1, X( I, 1 ), LDX )
  256.             CALL ZGEMV( 'Conjugate transpose', I-1, N-I+1, -ONE,
  257.      $                  A( 1, I ), LDA, X( I, 1 ), LDX, ONE, A( I, I ),
  258.      $                  LDA )
  259.             CALL ZLACGV( I-1, X( I, 1 ), LDX )
  260. *
  261. *           Generate reflection P(i) to annihilate A(i,i+1:n)
  262. *
  263.             ALPHA = A( I, I )
  264.             CALL ZLARFG( N-I+1, ALPHA, A( I, MIN( I+1, N ) ), LDA,
  265.      $                   TAUP( I ) )
  266.             D( I ) = ALPHA
  267.             IF( I.LT.M ) THEN
  268.                A( I, I ) = ONE
  269. *
  270. *              Compute X(i+1:m,i)
  271. *
  272.                CALL ZGEMV( 'No transpose', M-I, N-I+1, ONE, A( I+1, I ),
  273.      $                     LDA, A( I, I ), LDA, ZERO, X( I+1, I ), 1 )
  274.                CALL ZGEMV( 'Conjugate transpose', N-I+1, I-1, ONE,
  275.      $                     Y( I, 1 ), LDY, A( I, I ), LDA, ZERO,
  276.      $                     X( 1, I ), 1 )
  277.                CALL ZGEMV( 'No transpose', M-I, I-1, -ONE, A( I+1, 1 ),
  278.      $                     LDA, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
  279.                CALL ZGEMV( 'No transpose', I-1, N-I+1, ONE, A( 1, I ),
  280.      $                     LDA, A( I, I ), LDA, ZERO, X( 1, I ), 1 )
  281.                CALL ZGEMV( 'No transpose', M-I, I-1, -ONE, X( I+1, 1 ),
  282.      $                     LDX, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
  283.                CALL ZSCAL( M-I, TAUP( I ), X( I+1, I ), 1 )
  284.                CALL ZLACGV( N-I+1, A( I, I ), LDA )
  285. *
  286. *              Update A(i+1:m,i)
  287. *
  288.                CALL ZLACGV( I-1, Y( I, 1 ), LDY )
  289.                CALL ZGEMV( 'No transpose', M-I, I-1, -ONE, A( I+1, 1 ),
  290.      $                     LDA, Y( I, 1 ), LDY, ONE, A( I+1, I ), 1 )
  291.                CALL ZLACGV( I-1, Y( I, 1 ), LDY )
  292.                CALL ZGEMV( 'No transpose', M-I, I, -ONE, X( I+1, 1 ),
  293.      $                     LDX, A( 1, I ), 1, ONE, A( I+1, I ), 1 )
  294. *
  295. *              Generate reflection Q(i) to annihilate A(i+2:m,i)
  296. *
  297.                ALPHA = A( I+1, I )
  298.                CALL ZLARFG( M-I, ALPHA, A( MIN( I+2, M ), I ), 1,
  299.      $                      TAUQ( I ) )
  300.                E( I ) = ALPHA
  301.                A( I+1, I ) = ONE
  302. *
  303. *              Compute Y(i+1:n,i)
  304. *
  305.                CALL ZGEMV( 'Conjugate transpose', M-I, N-I, ONE,
  306.      $                     A( I+1, I+1 ), LDA, A( I+1, I ), 1, ZERO,
  307.      $                     Y( I+1, I ), 1 )
  308.                CALL ZGEMV( 'Conjugate transpose', M-I, I-1, ONE,
  309.      $                     A( I+1, 1 ), LDA, A( I+1, I ), 1, ZERO,
  310.      $                     Y( 1, I ), 1 )
  311.                CALL ZGEMV( 'No transpose', N-I, I-1, -ONE, Y( I+1, 1 ),
  312.      $                     LDY, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
  313.                CALL ZGEMV( 'Conjugate transpose', M-I, I, ONE,
  314.      $                     X( I+1, 1 ), LDX, A( I+1, I ), 1, ZERO,
  315.      $                     Y( 1, I ), 1 )
  316.                CALL ZGEMV( 'Conjugate transpose', I, N-I, -ONE,
  317.      $                     A( 1, I+1 ), LDA, Y( 1, I ), 1, ONE,
  318.      $                     Y( I+1, I ), 1 )
  319.                CALL ZSCAL( N-I, TAUQ( I ), Y( I+1, I ), 1 )
  320.             ELSE
  321.                CALL ZLACGV( N-I+1, A( I, I ), LDA )
  322.             END IF
  323.    20    CONTINUE
  324.       END IF
  325.       RETURN
  326. *
  327. *     End of ZLABRD
  328. *
  329.       END
  330.